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Abstract. Numerical MHD codes have become extraordinarily powerful tools 
with which to study accretion turbulence. They have been used primarily to 
extract values for the classical a parameter, and to follow complex evolutionary 
development. Energy transport, which is at the heart of classical disk theory, 
has yet to be explored in any detail. Further topics that should be explored 
by simulation include nonideal MHD, radiation physics, and outburst behavior 
related to the temperature sensitivity of the resistivity. 



1. Introduction: The Need for MHD Simulations 

Magnetic fields, particularly weak magnetic fields, strongly destabilize differen- 
tial rotational in astrophysical accretion disks, allowing them to live up to their 
name. The magnetorotational instability (MRI) is but one manifestation of a 
very general behavior of weak magnetic fields in plasmas: by forcing angular 
momentum and entropy transport to follow field lines, they allow free energy 
gradients (temperature and angular velocity) to become sources of instability. 
Despite the fact that they have only been recently understood in detail, the insta- 
bilities are not difficult to simulate and to study numerically. (A single processor 
workstation will do.) For this reason, numerical investigations of spontaneously 
arising, undriven accretion disk turbulence are based on the MRI. 

Nevertheless, the notion that simple differential rotation is intrinsically un- 
stable via nonlinear hydrodynamic instability (Dubrulle 1993) continues to at- 
tract adherents (e.g., Richard & Zahn 1999; Richard 2003), in part because 
there is as yet no mathematical proof of nonlinear stability. It is argued that 
shear layers are known to be linearly stable, but extremely sensitive to nonlin- 
ear disruption, thus Keplerian differential should also be nonlinearly unstable, 
provided that the Reynolds number is sufficiently high. 

In fact, one does not even have to resort to nonlinearity to find an example 
of shear induced instability. Linear global hydro instabilities can be present un- 
der some circumstances, provided that proper boundary conditions are specified. 
These "Papaloizou-Pringle" (Papaloizou & Pringle 1984; Goldreich, Goodman, 
&: Narayan 1986) instabilities, which involve communication by trapped waves 
on either side of a corotation radius, can disrupt systems that would ordinarily 
be stable by the classical Rayleigh criterion of increasing outward specific an- 
gular momentum. Numerical study of the PP instability poses no unresolvable 
difficulties, and it may be followed well into its nonlinear regime (Hawley 1991). 
Indeed, the emergent characteristic behavior of this instability is an instructive 
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example of how a controlled simulation can be used in pursuing an answer to 
a well-posed problem. The PP instability creates spiral structure in the disk, 
rather than turbulence, since this turn out to be a more efficient angular mo- 
mentum transport mechanism in a hydrodynamical flow. But both because the 
instability is very sensitive to the proper boundary conditions, and because the 
instability weakens dramatically when the pressure gradients vanish, the rele- 
vance PP instabilities to Kcplcrian disks is far from clear. 

More to the point, numerical simulations offer no support for the notion that 
local differential rotation is itself nonlinearly unstable. Proponents of nonlinear 
instability argue that this is simply because existing codes do not have a high 
enough effective Reynolds number, and are thus too diffusive. I must confess 
to being somewhat baffled by this argument, since the nonlinear instability of a 
shear layer is readily demonstrated, and it certainly does not require outlandishly 
high Reynolds numbers. Moreover, tight convergence has been achieved between 
codes with very different diffusive properties, and these show unambiguous non- 
linear stability of local Keplerian flow (Hawley, Balbus, Sz Winters 1999). 

It has been recently noted that the local "shearing box" equations of motion 
(e.g. Balbus &; Hawley 1998), which are fully nonlinear but do not include radial 
curvature terms, exhibit an exact scale invariance (Balbus 2004b). Specifically, 
if w(r, t), P{r,t), and p{r,t) represent exact solutions for the velocity, pressure 
and density respectively, then 

-^w(/3r,t), piPr,t), -^P{pr,t) 

where /3 is an arbitrary scaling factor, are also exact solutions. This means that 
Keplerian flow is no more unstable at very small scales then at larger, com- 
putationally accessible scales. Any putative local disk instability in Keplerian 
flow would have to exhibit behavior very different from that seen in the classical 
breakdown of planar shear flow. Clearly, Coriolis forces are strongly stabilizing, 
rigorously so in the linear regime, and it is not surprising that the stabiliza- 
tion extends into the nonlinear regime as well. As we have already noted, this 
does not mean that any rotationally-based instability is proscribed, it simply 
means that local shear is insufficient to disrupt Keplerian flow. Nonlinear in- 
stability buffs arc forced to argue that conditions beyond the reach of modern 
supercomputers arc required to trigger instability, a position that may become 
increasingly untenable as time goes on. As a very practical matter, the MRI 
remains the only instability on the market that manages to sustain long term, 
unforced turbulence in numerical simulations of Keplerian disks. In the end, 
numericists really don't have a choice. 

2. Tools of the Trade: Correlated Fluctuations 

Accretion disk turbulence cannot yet be directly observed in nature, and in the 

laboratory the MRI has only very recently been identified (Sisan et al. 2004). 
Three-dimensional local MHD simulations are almost ten years old at the time 
of this writing (Hawley, Gammie, 8z Balbus 1995). Three-dimensional global 
simulations followed shortly thereafter (Armitage 1998; Hawley 2000, 2001). 
In fact, two-dimensional global simulations date back two decades (Uchida & 
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Shibata 1985), well before MRI calculations appeared in local guise, and before 
the MRI had been identified in accretion disks. By now, global MHD simulations 
are a small industry, with several groups around the world investigating many 
different specialized problems. What can be learned from these efforts? 

We must first understand what it is that we don't understand, and for that 
we should begin with the tenets of classical disk theory (CDT). CDT is more 
widely known as viscous disk theory, a moniker that surely needs to be laid to 
rest. CDT helped to orient our intuitions during the early stages of disk theory, 
but now that this phase has past, it is time for disk viscosity to be jettisoned. 
The analogy has became more of a hindrance than anything else — in some 
papers enhanced viscosity and turbulence are viewed as competing processes! 
It is time for disk theorists to move out of the 1970s, and relegate "anomalous 
viscosity" to the cultural dustbin along with afros, leisure suits, and disco. 

2.1. Stress Tensor 

The fundamental transport quantity that is at the center of any theory of tur- 
bulent disk transport is the stress tensor. Unlike a viscous stress, the true stress 
tensor is entirely a property of the flow, not of the constituent fluid. It is defined 
in terms of correlated velocity fluctuations and correlated magnetic flclds. Let 
u be the difference between the local disk velocity and pure Kcplcrian rotation 
on cylinders. Let ua be the Alfven velocity associated with the magnetic field 
B, and mass density p, 

and is not necessarily a fiuctuating quantity. We assume that the magnitudes 
of Ua and u are comparable, and that both are much less that the unperturbed 
rotation velocity, RQ, where R is the cylindrical radius. The density and pressure 
fluctuations, 6p and 6P are assumed to be small compared with their mean 
values, p and P, respectively. The stress tensor has, in principle, six independent 
components, but in practice only the azimuthal components play an essential 
role in disk theory, 

T = p{{uu^ - UAUAcp)) (2) 

and we shall henceforth refer to these components as the stress tensor. Extract- 
ing an accurate mean value of T from numerical simulations requires patience, 
as very long time averaging may be required (Stcinacker & Papaloizou 2002). 

The stress tensor plays two mathematically related, but physically distinct 
roles is disk theory: (1) it is directly proportional to the outward angular mo- 
mentum flux through the fluid, and hence to the classical a parameter; and 
(2) it is directly proportional to the rate at which energy is extracted from the 
background differential rotation of the disk. The latter is the free energy source 
that powers the turbulent fluctuations themselves. The central tenet of CDT is 
that this energy goes nowhere. Instead, it is all locally dissipated and radiated 
away on the spot. The reader will recall the well-known factor of 3 discrepancy 
between the radiated energy and the local orbital energy that is released. This 
is usually explained as a nonlocal consequence of the presence of an energy flux, 
but I prefer to think of it as the result of the local injection of the free energy 
of differential rotation. In any case, this assumption of local energy dissipation 
is very powerful, and needs to be carefully checked. 
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2.2. Fluctuation Equations 

Hydrostatic equilibrium The steady state disk equations iiave been worked out 
under assumptions of section 2.1. The fundamental equation of hydrostatic 
equilibrium is 

m'^en = ivP + V$ (3) 

where ep, is a unit vector in the radial direction, and ^ is the central potential. 
Notice that hydrostatic equilibrium may be expressed in terms of mean values, 
and is not directly affected by the turbulent fluctuations. Equation and those 
to follow, are two-dimensional: the process of extracting the mean averages out 
the azimuthal structure. Hydrostatic equilibrium breaks down, of course, in 
regions of thermal or transonic flow. 

Mass conservation By way of contrast, transport processes are intimately 
linked to correlations in fluctuating quantities, and to the magnetic field. (This 
is true even if the field is weak.) The mass fiux is 

(pv) = pv+(5p(5v), (4) 

and in steady state 

V • (pv) = (5) 

The final fluctuating term in equation @ is generally comparable to the mean 
flow term, though it is dropped in viscosity analogues. 

Angular momentum Angular momentum conservation is expressed by the steady 
state equation 

V • {R^n{pv) + RT) = (6) 

The first term represents that angular momentum transported directly by the 
mass fiux, whereas the second is fiux passing through the disk gas. In CDT, 
only radial transport is considered, and T is replaced by a viscous stress. 

Energy conservation Energy is transported through the body of the disk via 
the thermal energy flux, 

Fth = ^(5v<5r) (7) 
7-1 

where 7 is the adiabatic index, and r is the normalized temperature kT/^ {k is 
the Boltzmann constant and /x the mean mass per particle). This form of Fth is 
very general. The velocity-temperature correlation is responsible for convective 
transport, as well as the energy transported by ordinary sound waves. The 
energy equation for a turbulent disk is (Balbus 2004a): 

- Tm^^ = V • Fth + r(pv) • VS + Qrad (8) 

where Tpi^ is the radial component of T, Qrad is the volumetric radiative loss 
rate, and S is in essence the mean entropy, 
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Equation Q states that the rate at which energy is locally extracted from the 
differential rotation is equal to the rate at which it is carried away by waves, 
dissipated as heat, and radiated away. 

In CDT, all terms on the right side of the equation are ignored except for 
radiative losses, Qrad- This is justified in the following sense. It is customary to 
express the stress tensor in terms of a simple dimensionless scaling, 



'-R(p 



aP (9) 



where a dimensionless parameter that embodies both the magnitudes of the 
velocity fluctuations as well as their relative degree of correlation. This will of 
course be immediately recognized as the relationship that gives "a-disk" theory 
its name. Since there is more than one correlation tensor of interest here, let us 
relabel the above alpha as ai, since T/j^ is associated with angular momentum 
transport. We are also free to invoke a scaling law of the form 

Fe = oePt, (10) 

where Fe is a radial energy flux. The neglect of this thermal energy flux will 
then be justified if 

The first ratio is an inverse Mach number for the rotation velocity, and will be 
small (by definition) for a thin Keplerian disk. This presumably includes CV 
disks. The energy flux term will then be negligible if q^; does not much exceed 
ol. While this is not a priori unreasonable, neither is it very well tested. It is 
the kind of question that should be asked of numerical simulations. 

2.3. Global Simulations 

Most of the the global disk simulations to date follow the following prototype. 
The gas is a polytropic torus, in the gravitational field of a "black hole." This 
means that the potential field is given by (Paczyhski & Wiita 1980): 

^ GM 

^ = -1 r, 12 

|r- rgl 

where Vg is the usual Schwarzschild (or gravitational) radius, 2GM/c^ . (Symbols 
have their usual meanings.) This initial state is chosen to be a torus because it 
allows a simple hydrostatic equilibrium to be constructed, but it is also physically 
reasonable: a reservoir of angular momentum bearing material located at a large 
distance from the origin. The equilibrium is of course nonmagnetic; as soon as 
a small field is added, the MRI quickly dominates the evolution. 

The original motivation for introducing the Paczyhski- Wiita potential was 
that it captured some key features of true black hole dynamics in a user-friendly 
Newtonian format. For the numericist, this is a particularly convenient potential 
function, since the inner boundary is pure outflow, or, from the point of view 
of the hole, an inflow at r = r^. The use of this potential when the central 
object is not a black hole is obviously a cheat, but it does flnesse the problem of 
how to handle the numerically intractable disk-star boundary layer. The price 
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to be paid is that all information of this observationally critical region is lost. 
But if the scope of the simulation is restricted to the study of the fundamental 
turbulent transport properties in the body of the disk, this approach makes some 
sense, particularly at this early stage. One must be willing to accommodate the 
possibility, however, that the existence of the inner boundary layer could cause 
global changes in ways that are difficult to foresee. 

If we take at face value the notion that numerical codes calculate the evo- 
lution of a strictly polytropic disk, equation ((HJ leaves only one choice for its 
dominant balance: 

This is a balance between the thermal energy flux divergence and the rate at 
which energy is exchanged with differential rotation, a purely adiabatic process. 
In the case of WKB waves, it leads to the conservation of wave action, rather 
than wave energy. More generally, this balance would involve extraction of the 
free energy of differential rotation followed by its active transport through the 
disk. This would involve an adiabatic coupling between evanescent and wavelike 
modes. Does this really happen in simulations, let alone real disks? 

It hasn't been checked. ZEUS-like codes, which base their energetics on an 
internal, as opposed to a total energy equation, tend to be lossy when comput- 
ing small scale flow structure. This is not necessarily bad, and in fact it mimics 
the behavior of real disks. In numerical tests of the stability of hydrodynamical 
flow, these energy losses were tracked and their rapid growth used as a hallmark 
for the presence of turbulence (Balbus, Hawley, &: Stone 1996). These losses in- 
troduce an effective Qrad term in the calculation, even though it is not explicitly 
included in the code. The interesting question is how much free energy is locally 
"radiated" and how much is transported elsewhere, since this will determine 
whether local disk models make sense. The fact that simulations do not lead 
to very thin disks suggests that local grid losses are not overwhelming, but a 
detailed accounting should be part of the numerical diagnostics. The existence 
of thin disks as well as the local R~^^^ scaling law seen in some eclipsing bina- 
ries like Z Cha (Frank, King, & Raine 2002) are certainly consistent with local 
dissipation, but it is not obvious that they require it. 



3. Disk Morphology 

Disk morphology predictions may be more robust than is often recognized. The 
angular velocity Q tends to be constant on cylinders, since z-dependent rotation 
profiles are unstable (Goldreich & Schubert 1967). If so, the right hand side of 
equation (jSJ is a pure gradient, and whatever the actual constitutive relationship 
between P and p may be, the equilibrium functional relationship between these 
two quantities will be that isobaric and isochoric surfaces coincide. Moreover, 
these surfaces are also equipotentials, provided that the effective centrifugal 
potential is included (e.g., Frank et al. 2002). 

In principle, the rotation profile ^1{R) is a free functional parameter, but 
in practice almost all MRI simulations show that a Keplerian power law rapidly 
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emerges: 



GM cos/3 

^3 



(14) 



where cos/9 is a proportionality constant not very different from unity. The 
reason for this seems to be connected with the vigorous outward angular mo- 
mentum transport that always accompanies the MRI in disks, which spreads 
the initial torus radially, diluting the dynamical effects of pressure gradients in 
the process. The disk then seems to follow a simple a quasi-Keplerian profile, 
linear in GM, with residual pressure support allowing rotation at slightly be- 
low its pressure-free value. For a polytropic equation of state, P = Kp^ , the 
equilibrium density profile satisfies 



where Woo is the enthalpy at infinity, and A is the colatitude, 7r/2 — 9, where 9 
is the usual spherical angle. An isothermal equation of state yields: 



There is a qualitative difference between an adiabatic and an isothermal 
disk. An adiabatic disk has a well-defined sharp edge at the value of A corre- 
sponding to p = 0, whereas an isothermal disk trails off exponentially above the 
disk boundary A = /3. The cquipotential surfaces themselves satisfy the equation 



Contours for a fiducial example are shown in Figure 1. When C > 0, the 
contours are open, extending to infinity. When C < 0, they curl back through 
the midplane and reconnect at the origin. The critical C = contour is the 
straight line X = P, dividing the two classes. Notice that all contours converge 
to this critical contour as r ^ 0. 

The zone of "contour convergence" seems to be the site of the launching of 
a jet-like outflow in numerical simulations, giving some credence to this analytic 
model. (The jet appears to be thermal, rather than magneto-centrifugal, in 
origin.) The tight grouping of the equipotential surfaces in the convergence 
zone allows dissipative processes to move fluid elements from one equipotential 
surface to another fairly easily in this region. Incoming disk material cannot be 
swallowed as fast as it arrives, and is forced onto upward onto open equipotential 
surfaces. The inner most allowable surface will always coincide more closely to 
a constant angular momentum cylinder then will surface closer to the disk. This 
will be the easiest contour for the fluid to flow along as it leaves the disk, requiring 
a minimum of angular momentum change. And indeed, the simulations show a 
pile-up of material close to the last open equipotential surface, confined against 
this surface from the inside by a low density magnetic corona. Significantly, the 
isothermal runs do not show the formation of a jet. This can be understood both 
because the softer equation of state allows material to be more easily swallowed 
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isoenergy (density, pressure) plots for H >0, (3=5i/4 




X H /(G M) 

oo 

Figure 1. Equipotential adiabatic disk contours for the case when is 0.84 
of its Keplerian value, corresponding to an opening wedge angle of 45° , marked 
by the separatrix line marked OUTFLOW. Open contours become very closely 
packed near the p = boundary. 

at the center of flow, and because the density is exponentially curtailed above 
the cone X = (3. 

Figure 2 shows a detailed simulation by De Villiers, Hawley, & Krolik 
(2003). The plot displays azimuthally averaged density contours after a pe- 
riod of time corresponding to 10 orbital periods at the pressure maximum of the 
initial torus used in the simulation. This is a general relativistic calculation of 
an adiabatic gas, but except very near the "plunging region", the gravitational 
potential is nearly Newtonian. The tear drop cross section predicted by analytic 
theory and inner jets are both present. 

4. Beyond Ideal MHD 

4.1. Hall and Ohmic Processes 

At temperatures below 2000 K, the level of ionization in a CV disks can fall 
below the minimum needed to keep the gas in the ideal MHD regime. Assuming 
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Figure 2. Cross section of density contours of an adiabatic accretion disk 
surrounding a Kerr black hole. Contour shape and presence of jet-like outflow 
follow equipotential surfaces from equation (17). 



that dust grains are unimportant, the relative mutual drift velocities of ions, 
electron, and neutrals gives rise to ambipolar diffusion, Hall electromotive forces 
(HEMFs), and Ohmic dissipation. Ambipolar diffusion is important in interstel- 
lar clouds, but for the much denser disks of interest here, HEMFs and Ohmic 
losses are of primary interest. Numerical MHD has recently begun to study 
these nonideal effects, which can have a decisive influence over disk behavior. 

Ohmic losses were studied in local disk calculations by Fleming & Stone 
(2003), who found that MHD turbulence could be maintained only if the mag- 
netic Reynolds number 

ReM = -p:, 

7]il 

where r/ is the resistivity, exceeded a critical value, ReM{crit). The value of 
Reuicrit) depends upon the field geometry. For runs with a mean vertical field 
present, Fleming & Stone found critical values of about 100; these increased to 
10'* for runs without a mean vertical field. Once Rcm < ReMicrit), turbulence 
is quickly suppressed. The addition of HEMFs (Sano & Stone 2002) seems not 
to alter very much the values of ReM{crit), but the saturations levels are more 
sensitive. These can be significantly raised by the presence of the Hall effect. 

Numerical non-ideal MHD offers yet another venue of opportunity. The 
ohmic resistivity is extraordinarily temperature sensitive in the transition regime. 
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because it depends upon the electron ionization fraction. If thermal in origin, 
this fraction depends upon Boltmann factors of the form exp(— //feT) where 
the ionization energy / is well above kT, which accounts for the temperature 
sensitivity. Small temperature fluctuations cause large changes in the ionization 
fraction. An increase in T leads to a decrease in rj, and, it may be expected, to 
higher levels of turbulence. This is a prescription for yet higher temperatures, 
and a possible runaway. A similar instability is associated with T decreases and 
suppression of turbulence. How this resolves itself is not yet known, but it seems 
to mc that the problem is within the computational realm. 

Applications of this behavior to DN eruptions have been explored (Gammie 
& Menou 1998; Menou 2000; Sano &: Stone 2003), but interpreting the quiescent 
stage in this context remains problematic. Fleming &; Stone (2003) make the 
interesting and provocative suggestion that magnetically "dead" zones may be 
affected by surrounding lower density magnetically active regions. In this view, 
the dead zones may be only morbid, and host an active stress tensor maintained 
by outwardly transporting trailing density waves. If this is in fact the way that 
such low ionization regions behave, it is a bit of a coup for numerical studies of 
accretion disks. It is difficult to imagine a compelling a priori analytic argument, 
for example, that would give life to the dead zone in this way. 

The difficulty with the low state has been to avoid having the a parameter 
drop all the way to zero. Gammie &: Menou (1998) suggested global hydro- 
instabilities, Menou (2000) tidal interactions, and Sano & Stone (2003), in their 
numerical calculation, treated the low state "by hand". The view put forth 
above offers another possibility: the quiescent disk state might be a hybrid wave- 
turbulent system, while the high state is a flow fully governed by the MRI. In its 
low state, the disk is obstructed by a sluggish but extended central zone, while 
accretion occurs more easily in a rarified atmosphere around the low ionization 
region. Even in the low hybrid state, the primary accretion process may still be 
largely MHD, but affecting only the tenuous gas. 

4.2. Radiation 

Radiation is essential to CDT, since this is the fate of all the free energy that 
is locally dissipated by turbulence. Very little work has been done, however, 
to include radiative losses in numerical simulations. What has been done, has 
focused on the case where the radiation energy rivals or dominates the thermal 
energy density, which can happen in the inner regions of black hole accretion. 
The linear stability of a magnetized, stratified, radiative gas has been investi- 
gated by Blaes &; Socrates (2001). Despite the complexity of the full problem, 
the MRI emerges at the end of the day unscathed, its classical stability criterion 
dU^/dR > remaining intact. 

The proper interpretation of a linear analysis is slightly unclear, since the 
unperturbed state is already presumably fully turbulent due to the MRI; a ro- 
tationally stable disk lacks an internal energy source. The interplay between 
MHD turbulence and radiation is intrinsically nonlinear. Turner, Stone, & Sano 
(2002) have studied a local, radiative, axisymmetric, shearing box. The linear 
calculations of Blaes Sz Socrates were confirmed in detail, and the nonlinear flow 
fully developed. As in standard MRI simulations, the stress is dominated by the 
Maxwell component, which is a factor of a few larger than the Reynolds terms. 
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Clcarlv; much could be learned by returning to the case where radiation is 
an important loss mechanism, but does not dominate the local energy density — 
i.e., the Shakura & Sunyaev (1973) prototype. The technical expertise brought 
to bear on radiation dominated disks would serve this simpler problem well. Es- 
tablishing that a radiative shearing box does, in fact, evolve toward something 
resembling the classical Shakura-Sunyaev paradigm would be an important mile- 
stone for accretion theory. 



5. Summary 

• The MRI remains the leading candidate for the origin of enhanced angular 
momentum and energy transport in accretion disks. A compelling case for 
sustained hydro dynamical turbulence has yet to be made. 

• Despite very real limitations, the most serious of which is probably still 
the limited dynamical range, numerical simulations have become truly 
powerful. They have allowed us to follow the global evolution of polytropic 
disks, and to witness the build-up of a Keplerian disk, corona, and central 
jet — starting from nothing but a simple constant angular momentum torus. 
Many of these features can be understood analytically, at least in a crude 
sense. 

• Energetics (including radiation and dissipative physics) remains in its 
infancy-we know less than we think. Numerical disk models have yet 
to be placed on the observational plane. 

• We have not tested under what conditions the fundamental formula of 
phenomenological disk theory 

Q- = -Tr4> 



dlnR 



is valid. Even the eclipse mapping observations that show an power 
law dependence in the temperature profile of CV disks do not require 
this assumption to be valid, though they are certainly consistent with 
it. Investigating this question numerically should be a central goal for 
numericists. 

A closely related point: We arc completely ignorant of the role of ther- 
mal energy transport, 6v6T. It need not be small relative to rotational 
transport, and can be extracted from simulations. 

The consequences for the MRI of departures from ideal MHD in the form 
of Hall EMFs and ohmic resistance are critical to understand. Only one 

team, Sano & Stone, has studied this mimcrically, and then only on fairly 
coarse grids. Resolution may be key here. The question that needs to be 
thoroughly explored is under what circumstances does the MRI turn of? 

The temperature dependence of the resistivity has not been touched nu- 
merically. Clearly, there is a real potential to explore eruptive behavior 
here. 
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• For all their foibles and limitations, numerical MRI simulations have taken 
their place at the helm of theoretical accretion disk studies, and they are 
likely to remain there for the foreseeable future. 
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